
\chapter{线性方程组}
线性方程组是一个非常常见的问题。
很多看似不相干的问题也可以转化为线性方程组求解问题。
比如常微分方程：
\begin{equation}
\frac{d^2}{dx^2} u(x) + p(x) \frac{d}{dx} u(x) +  q(x) u(x) = r(x),
\label{derivative-equation-order2}
\end{equation}
在$[a,b]$上成立，另有边界条件$u(a) = u(b) = 0$。

如果在$[a,b]$上取$n$个节点，将各阶导数用有限差分公式代替，
\begin{eqnarray}
\frac{d^2}{dx^2} u(x) |_{x_i} &\approx& ( f(x_{i+1}) + f(x_{i-1}) - 2f(x_i) )/h^2,
\label{double-derivative-finite-difference}
 \\
\frac{d}{dx} u(x)|_{x_i} &\approx& ( f(x_{i+1}) - f(x_{i-1}) )/(2h)
\label{derivative-finite-difference}
\end{eqnarray}
带入微分方程，就会得到$n$个线性方程，变量为$u(x_1), u(x_2), \cdots, u(x_n)$，解出来，就得到了微分方程的一个数值解。

在这一章中，我们介绍常规的线性方程组解法：高斯消元法、LU分解法、三对角矩阵追赶法、迭代法。

\section{高斯消元法}
在小学三年级或者四年级，我们就学过了高斯消元法。
线性方程组如下
\begin{eqnarray}
a_{11} x_1 + a_{12} x_2 + \cdots + a_{1n} x_n &=& b_1 \nonumber\\
a_{21} x_2 + a_{22} x_2 + \cdots + a_{2n} x_n &=& b_2 \nonumber\\
\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots &=& \cdots \nonumber\\
a_{n1} x_1 + a_{n2} x_2 + \cdots + a_{nn} x_n &=& b_n 
\end{eqnarray}
或者记作
\begin{equation}
A \vec{x} = \vec{b},
\end{equation}
其中
\begin{eqnarray}
A = \begin{pmatrix}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots & \vdots & \vdots & \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn}
\end{pmatrix},
~~~~
\vec{x} = \begin{pmatrix}
x_1\\
x_2\\
\vdots\\
x_n
\end{pmatrix},
~~~~
\vec{b} = \begin{pmatrix}
b_1\\
b_2\\
\vdots\\
b_n
\end{pmatrix}
\end{eqnarray}
这个方程有唯一解的前提是，$A$是满秩矩阵，即$|A| \neq 0$。

我们注意到，如果交换A的两行，并相应交换$\vec{b}$，不影响解$\vec{x}$；如果将A的任一行乘以一个系数再加到另一行上，并对$\vec{b}$做相应操作，不影响解$\vec{x}$。
如果交换A的两列，得到的解则与$\vec{x}$不一样，相当于$\vec{x}$交换相应的两行。

我们可以将A与$\vec{b}$合在一块，记作$\mathcal{A}$，
\begin{equation}
\mathcal{A} \equiv \begin{pmatrix}
A & \vec{b}
\end{pmatrix}
\end{equation}
线性方程组则可以表达为
\begin{eqnarray}
\mathcal{A} \begin{pmatrix}
\vec{x} \\ -1 \end{pmatrix} = 0
\end{eqnarray}

对于$\mathcal{A}$中的第$j=2,\cdots,n$行，我们做行变换，
\begin{eqnarray}
\mathcal{A}_{jk} = - \frac{ \mathcal{A}_{j1}}{ \mathcal{A}_{11}} \mathcal{A}_{1k} + \mathcal{A}_{jk}, ~~k=1,2,\cdots, n,
\label{Gauss1}
\end{eqnarray}
即用$- \frac{ \mathcal{A}_{jk}}{ \mathcal{A}_{11}}$乘以第一行，再加到第$j$行上。
这么做的后果是，第2-n行第1个元素等于零。
注意，上面的公式不是一个等式，而是一个算法操作，等号右边的$\mathcal{A}$是操作前的矩阵，等号左边的$\mathcal{A}$是更新后的矩阵。

类似地，我们可以将第3-n行第2个元素消为零，$\forall j = 2,3,\cdots,n$，
\begin{eqnarray}
\mathcal{A}_{jk} = - \frac{ \mathcal{A}_{j2}}{ \mathcal{A}_{22}} \mathcal{A}_{2k} + \mathcal{A}_{jk}, ~~k=2,\cdots, n,
\label{Gauss2}
\end{eqnarray}

如此这般，进行下去，第$i$次行变换为$\forall j = i,i+1,\cdots,n$，
\begin{eqnarray}
\mathcal{A}_{jk} = - \frac{ \mathcal{A}_{j,i}}{ \mathcal{A}_{i,i}} \mathcal{A}_{i,k} + \mathcal{A}_{j,k}, ~~k=i,\cdots, n,
\label{Gaussi}
\end{eqnarray}
最终$\mathcal{A}$会变成一个上三角矩阵：
\begin{equation}
\mathcal{A} = \begin{pmatrix}
\mathcal{A}_{11} & \mathcal{A}_{12} &  \cdots & \mathcal{A}_{1,n} & \mathcal{A}_{1,n+1} \\
0                & \mathcal{A}_{22} &  \cdots & \mathcal{A}_{2,n} & \mathcal{A}_{2,n+1} \\
0                & 0                &  \cdots & \mathcal{A}_{3,n} & \mathcal{A}_{3,n+1} \\
\vdots           & \vdots           &  \vdots & \vdots            & \vdots \\
0                & 0                &  \cdots & \mathcal{A}_{n,n} & \mathcal{A}_{n,n+1} 
\end{pmatrix}
\end{equation}
变换到这一步，求解 $\mathcal{A} \begin{pmatrix}
\vec{x} \\ -1
\end{pmatrix}
=0$ 就会变得非常简单：$\forall j = n, n-1, \cdots, 1$，
\begin{eqnarray}
x_j = \frac{1}{ \mathcal{A}_{j,j} } \{ \mathcal{A}_{j,n+1} - \sum^n_{k=j+1} \mathcal{A}_{jk} x_k \}
\end{eqnarray}

但是有一个地方可能出现问题，比如（\ref{Gauss1}）或（\ref{Gauss2}）中，如果$\mathcal{A}_{11}$或$\mathcal{A}_{22}$等于零，或者是一个非常小的双精度数，就无法操作，或者带来很大的误差。
所以我们需要在每次行变换之前，增加一步：
\begin{shaded*}
对第$k,\cdots,n$行进行行变换之前，如果$\mathcal{A}_{k,k}$的绝对值太小，就寻找$\mathcal{A}_{k,k}, \mathcal{A}_{k+1,k}, \cdots, \mathcal{A}_{n,k}$中的最大值，并交换相应的行与第$k$行。
\end{shaded*}
这就叫做寻找主元，相应地，这样做叫做高斯主元消去法。

如果我们把加减乘除看做是一次“操作”，那么如果不算换行操作，高斯消元法解一个n元线性方程组需要$O(n^3)$次操作。
所以，如果$n=1000$，电脑主频在$GHz$数量级，解这样的常微分方程组可能需要几秒钟。

下面的代码实现高斯主元消去法。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>

/*
 * GaussLinear: solves n-dim linear equations
 * The solution is stored in x
 * Permutation of rows is activated, when the principle element is smaller than min
 * If A is singular, a permutation will fail and function exit(1), signals error
 */
void GaussLinear(int n, double **A, double *x, double min){

        int i,j,k;
        double y;
        for(i=0;i<n;i++){
                if( fabs( A[i][i] ) < min ){// permutation is activated
                        for(j=i+1;j<n;j++){
                                if( fabs( A[j][i] ) > min ){
                                        for(k=i;k<n;k++){
                                                y = A[i][k];
                                                A[i][k] = A[j][k];
                                                A[j][k] = y;
                                        }
                                        break;
                                }
                        }
                        if(j==n){// didn't find any element larger than min
                                cerr<<" A is singular in GaussLinear(n,A,x,min).\n";
                                exit(1);
                        }
                }
                for(j=i+1;j<n;j++){// row operations
                        y = - A[j][i] / A[i][i];
                        for(k=0;k<n+1;k++){
                                A[j][k] = y * A[i][k] + A[j][k];
                        }
                }
                // for test only
                /*
                cout<<"i="<<i<<"\t A:"<<endl;
                for(j=0;j<n;j++){
                        for(k=0;k<n+1;k++){
                                cout<<A[j][k]<<",";
                        }
                        cout<<endl;
                }
                */
        }
        for(j=n-1;j>=0;j--){
                y=0;
                for(k=j+1;k<n;k++) y += A[j][k]*x[k];
                x[j] = 1/A[j][j] * ( A[j][n] - y );
        }
}

int main(){

        int i,j,n=3;
        double B[3][4] = {{1,1,1,3},{1,2,2,5},{1,3,4,8}};
        double **A = new double * [n];
        for(i=0;i<n;i++){
                A[i] = new double [n+1];
                for(j=0;j<n+1;j++){
                        A[i][j] = B[i][j];
                }
        }
        double *x = new double [n];
        GaussLinear(3,A,x,1E-4);
        cout<<"x=";
        for(int i=0;i<3;i++) cout<<x[i]<<",";
        cout<<endl;
        delete [] x;
        for(i=0;i<n;i++) delete [] A[i];
        delete [] A;
        return 0;
}
\end{lstlisting}

\section{LU分解}
下三角矩阵(L矩阵)就是这样的：
\begin{equation}
L = 
\begin{pmatrix}
\begin{matrix}
\alpha_{00} & \\
\alpha_{10} & \alpha_{11}
\end{matrix}
& \text{\Large 0} \\
\begin{matrix}
\vdots & \vdots \\
\alpha_{n-1,0} & \alpha_{n-1,1}
\end{matrix}
&
\begin{matrix}
\ddots & \\
\cdots & \alpha_{n-1,n-1}
\end{matrix}
\end{pmatrix}
\end{equation}

而上三角矩阵（U矩阵）是这样的：
\begin{equation}
U = 
\begin{pmatrix}
\begin{matrix}
\beta_{00} & \beta_{01} \\
 & \beta_{11}
\end{matrix}
& \begin{matrix}
\cdots & \beta_{0,n-1} \\
\cdots & \beta_{1,n-1}
\end{matrix} \\
\text{\Large 0} 
&
\begin{matrix}
\ddots & \vdots \\
 & \beta_{n-1,n-1}
\end{matrix}
\end{pmatrix}
\end{equation}

任意$n\times n$的矩阵A可以分解为L矩阵与U矩阵的乘积:
\begin{eqnarray}
A &=&
\begin{pmatrix}
\begin{matrix}
a_{00} & a_{01} \\
a_{10} & a_{11}
\end{matrix}
& \begin{matrix}
\cdots & a_{0,n-1} \\
\cdots & a_{1,n-1}
\end{matrix} \\
\begin{matrix}
\vdots & \vdots \\
a_{n-1,0} & a_{n-1,1}
\end{matrix}
&
\begin{matrix}
\ddots & \vdots \\
\cdots & \alpha_{n-1,n-1}
\end{matrix}
\end{pmatrix}
\nonumber\\
&=& \begin{pmatrix}
\begin{matrix}
\alpha_{00} & \\
\alpha_{10} & \alpha_{11}
\end{matrix}
& \text{\Large 0} \\
\begin{matrix}
\vdots & \vdots \\
\alpha_{n-1,0} & \alpha_{n-1,1}
\end{matrix}
&
\begin{matrix}
\ddots & \\
\cdots & \alpha_{n-1,n-1}
\end{matrix}
\end{pmatrix}
\begin{pmatrix}
\begin{matrix}
\beta_{00} & \beta_{01} \\
 & \beta_{11}
\end{matrix}
& \begin{matrix}
\cdots & \beta_{0,n-1} \\
\cdots & \beta_{1,n-1}
\end{matrix} \\
\text{\Large 0} 
&
\begin{matrix}
\ddots & \vdots \\
 & \beta_{n-1,n-1}
\end{matrix}
\end{pmatrix} \nonumber\\
\end{eqnarray}

{\bf Crout}算法为：
\begin{shaded*}
\begin{itemize}
\item[1] $\forall i=0,\cdots,n-1, ~~\alpha_{ii} = 1$。
\item[2] $\forall j=0,1,2,...,n-1$，即每一列，
        \begin{itemize}
        \item [i] $\forall i=0,\cdots,j, ~~\beta_{ij} = a_{ij} - \sum^{i-1}_{k=0} \alpha_{ik} \beta_{kj}$,
        \item [ii] $\forall i=j+1, \cdots, n-1,~~ \alpha_{ij} = \frac{1}{\beta_{jj}} ( a_{ij} - \sum^{j-1}_{k=0} \alpha_{ik} \beta_{kj} )$
        \end{itemize}
\end{itemize}
\end{shaded*}
像上面这样做，有以下特点：
\begin{itemize}
\item [1] 每一步所需要的$\alpha, \beta$值都已经被前面的步骤计算好了。
\item [2] 每一步产生的$\alpha, \beta$值都供后面的步骤使用。
\item [3] 每个$a_{ij}$都只被调用一次，在计算$\alpha_{ij}$或$\beta_{ij}$时被调用。
\item [4] 计算第$j$列时，如果调换$A,L,U$第$j$行与第$i$行（$i>j$），不影响$L,U$矩阵前面第$0,1,\cdots,j-1$列的结果，这些结果仍然是正确的。
\end{itemize}
根据前3个特点，在代码中只需要一个$n\times n$的矩阵A，
在计算中不断用新的值覆盖旧的值，最终，A的下三角部分（不包括对角元）即$L$矩阵（其对角元均为1），上三角部分即$U$矩阵。

如果直接这么做，有1个漏洞，如果算法2.i中得到的$\beta_{jj}$等于0或者非常小，2.ii可能引入很大的误差。
所以在2.ii中，需要加入行变换，确保$\beta_{jj}$足够大。
上面的特点4保证了，引入的行变换不影响已经完成的结果。
引入行变换以后，Crout算法如下：
\begin{shaded*}
\begin{itemize}
\item[1] $\forall i=0,\cdots,n-1$，计算$max(i) = |a_{ij}|_{max}, j=0,1,\cdots,n-1$。
\item[2] $\forall j=0,1,2,...,n-1$，即每一列，
        \begin{itemize}
        \item [i] $\forall i=0,\cdots,j, ~~a_{ij} = a_{ij} - \sum^{i-1}_{k=0} a_{ik} a_{kj}$,
        \item [ii] $\forall i=j+1, \cdots, n-1,~~ \alpha_{ij} =  a_{ij} - \sum^{j-1}_{k=0} a_{ik} a_{kj} $
        \item [iii] 比较\\
        $\{ |a_{jj}| / max(i), |a_{j+1,j}|/max(j+1), \cdots, |a_{n-1,j}|/max(n-1) \}$，\\
        若$i=imax$时$|\alpha_{ij}| / max(i)$最大，则交换第j行与imax行，并记下这个行变换。
        \item [iv] $\forall i=j+1, \cdots, n-1, a[i][j] /= a[j][j]$。
        \end{itemize}
\end{itemize}
\end{shaded*}
其中$max(i)$是每一行各元素绝对值的最大值，所以$\beta_{ij} / max(i)$是归一化以后的值。
做这个归一化，是因为有可能有一行所有数都非常小，我们关心的是相对的大小。

上面这些步骤做完以后，得到的并不是A的LU分解，而是A经过那些行交换（算法中已经逐一记下）以后的矩阵的LU分解。

做完LU分解以后，线性方程组求解就很容易了
\begin{equation}
A \vec{x} = L U \vec{x} = \vec{b},
\end{equation}
我们可以记
\begin{equation}
\vec{y} \equiv U \vec{x},
\end{equation}
则有
\begin{equation}
L \vec{y} = \vec{b},
\end{equation}
这很容易解，尤其是$L_{ii}=1,i=0,\cdots,n-1$，
\begin{shaded*}
$\forall j = 0,1,\cdots,n-1$，
\begin{equation}
y_j = b_j - \sum^{j-1}_{k=0} L_{jk} y_k,
\end{equation}
\end{shaded*}
解出$\vec{y}$以后，很容易根据$U\vec{x} = \vec{y}$得到$\vec{x}$，
\begin{shaded*}
$\forall j = n-1, n-2,\cdots,0$，
\begin{equation}
x_j = \frac{1}{ U_{jj} }( y_j - \sum^{n-1}_{k=j+1} U_{jk} y_k ),
\end{equation}
\end{shaded*}

在下面的代码中，ludcmp对系数矩阵A进行LU分解，并记下发生的行交换，ludcmp\_solve则用ludcmp的结果对$A\vec{x} = \vec{b}$进行求解。
这两个函数都参考自\cite{NumericalRecipes}，做了少量改动。

\begin{lstlisting}
/*
 * ludcmp():    LU decomposition of a non-singular matrix A
 * The L matrix has L_{ii} = 1, and other elements are written into lower-left part of a, while the U matrix is written into the upper-right part of a.
 * Note that ludcmp does not rule out all singular cases. Some singular A can also be successfully decomposed by ludcmp.
 * n:           dimension
 * a:           n x n matrix
 * pivot:       records pivotting
 * sign:        1/-1 for even/odd times of pivotting
 */
void ludcmp(int n, double **a, int *pivot, int & sign){

        int i,j,k,imax;
        double big, y;
        double *max = new double [n];
        for(i=0;i<n;i++){
                big = 0;
                for(j=0;j<n;j++){
                        if( fabs(a[i][j]) > big ){
                                big = fabs( a[i][j] );
                        }
                }
                if( big < 1E-9 ){
                        cout<<"a is singular in ludcmp(n,a,pivot,sign).\n";
                        exit(1);
                }
                max[i] = big;
        }
        //cout<<"max:"; for(i=0;i<n;i++)cout<<max[i]<<","; cout<<endl;

        sign =1;
        for(j=0;j<n;j++){
                for(i=0;i<=j;i++){
                        y=0;
                        for(k=0;k<= i-1;k++){
                                y+= a[i][k] * a[k][j];
                        }
                        a[i][j] -= y;
                }
                imax = j;
                big = fabs( a[j][j] ) / max[j];
                for(i=j+1;i<n;i++){
                        y=0;
                        for(k=0;k<=j-1;k++){
                                y+= a[i][k] * a[k][j];
                        }
                        a[i][j] -= y;
                        if( y = fabs( a[i][j] ) / max[i] > big ){
                                big = y;
                                imax = i;
                        }
                }
                if( imax != j ){
                        pivot[j] = imax;
                        sign *= -1;
                        y = max[j];// exchange max[j] and max[imax]
                        max[j] = max[imax];
                        max[imax] = y;
                        for(k=0;k<n;k++){// exchange line imax and line j
                                y = a[j][k];
                                a[j][k] = a[ imax ][k];
                                a[ imax ][k] = y;
                        }
                }
                else{
                        pivot[j] = j;
                }
                for(i=j+1;i<n;i++) a[i][j] /= a[j][j];
                //for test only
                /*
                cout<<"~~~~~~~~~ j="<<j<<endl;
                for(i=0;i<n;i++){
                        for(k=0;k<n;k++){
                                cout<<a[i][k]<<",";
                        }
                        cout<<endl;
                }
                */
        }
        delete [] max;
}

/*
 * ludcmp_solve:  solves linear equation Ax=b, the answer x is written into b
 * a is not A, a is LU decomposition of pivotted version of A, and the pivotting is recorded in pivot
 * When U is singular, throws out an error and exit the program.
 */
void ludcmp_solve(int n, double ** a, int * pivot, double * b){

        int j,k;
        double y;

        // do the pivotting on b
        for(j=0;j<n-1;j++){
                k = pivot[j];
                y = b[j];
                b[j] = b[k];
                b[k] = y;
        }
        //cout<<"b after pivoting: "; for(j=0;j<n;j++)cout<<b[j]<<","; cout<<endl;
        // calculate y = Ux, and write it onto b
        for(j=0;j<n;j++){
                y=b[j];
                for(k=0;k<=j-1;k++){
                        y -= a[j][k] * b[k];
                }
                b[j] = y;
        }
        //cout<<"y:"; for(j=0;j<n;j++)cout<<b[j]<<","; cout<<endl;
        // calculate x, and write it onto b
        for(j=n-1;j>=0;j--){
                y=b[j];
                for(k=j+1;k<n;k++){
                        y -= a[j][k] * b[k];
                }
                if( fabs(a[j][j])<1E-9 ){
                        cout<<"error: U is singular in the LU decomposition of A.\n";
                        exit(1);
                }
                b[j] = y/a[j][j];
        }
}
\end{lstlisting}

\paragraph{行列式与逆矩阵} 有了LU分解，行列式马上就得到了，因为
\begin{equation}
|A| = | LU | = |L| \cdot |U|,
\end{equation}
而$|L|=1$（因为其对角元都是1），并且
\begin{equation}
|U| = \sum^{n-1}_{i=0} U_{ii},
\end{equation}
所以得到
\begin{equation}
|A| = \sum^{n-1}_{i=0} U_{ii}.
\end{equation}
如果分解的是A经过行交换以后的矩阵，那么其行列式乘以$-1/1$（奇数次/偶数次行交换），就得到A的行列式。

逆矩阵则可以如下求得
\begin{equation}
A A^{-1} = LU A^{-1} = E,
\end{equation}
把$A^{-1}$分为$n$个列矩阵，E也分为$n$个列矩阵，进行线性方程组求解，就得到了$A^{-1}$。
当然，在这样做之前，需要判断A的行列式是否为零。

\section{三对角方程组}
三对角方程组为如下形式
\begin{equation}
\begin{pmatrix}
b_1 & c_1 &    &  & \\
a_2 & b_2 & c_2 &  & \\
    & \ddots & \ddots & \ddots \\
    &     & a_{n-1} & b_{n-1} & c_{n-1} \\
    &     &       & a_n     & b_n
\end{pmatrix}
\begin{pmatrix}
u_1 \\
u_2 \\
\vdots \\
u_{n-1} \\
u_n
\end{pmatrix}
= 
\begin{pmatrix}
r_1 \\
r_2 \\
\vdots \\
r_{n-1} \\
r_n
\end{pmatrix}
\end{equation}
比如，方程（\ref{derivative-equation-order2}）做（\ref{double-derivative-finite-difference}-\ref{derivative-finite-difference}）处理以后，得到的就是一个三对角方程组。

我们只讨论{\bf 对角占优}的三对角矩阵，这也是我们处理的问题经常拥有的属性，即
\begin{equation}
|b_j| > |a_j| + |c_j|, ~~~ j=0,\cdots,n-1
\end{equation}
前人证明了，在这个条件下，我们下面介绍的追赶法不会遇到零主元。

如果用高斯消元法或者LU分解法，当然可以解决上面的方程组，但未免有些杀鸡用牛刀的感觉。

高斯消元法与LU分解都需要$O(n^3)$次操作。

如果用“追赶法”，则只需要$O(n)$次操作。
\begin{shaded*}
\begin{itemize}
\item [1] 令 $u_1 = r_1/b_1$，$\gamma_1 = 0$，$\forall j = 2,3,\cdots,n$，做如下运算
            \begin{itemize}
            \item [] $\gamma_j = \frac{ c_{j-1} }{ b_{j-1} - a_{j-1} \gamma_{j-1} }$
            \item [] $u_j = \frac{ r_j - a_j u_{j-1} }{ b_j - a_j \gamma_j }$
            \end{itemize}
\item [2] $\forall j = n-1, \cdots, 1$,
            \begin{itemize}
            \item [] $u_j = u_j - \gamma_{j+1} u_{j+1}$
            \end{itemize}
\end{itemize}
\end{shaded*}
第1步向前迭代，第2步向后迭代。

下面我们证明最终结果是方程组的解。
第1步向前迭代以后，得到的$u_j$我们可以改记作为$u'_j$，则有$u'_1 = r_1/b_1$，
\begin{eqnarray}
\gamma_j &=& \frac{ c_{j-1} }{ b_{j-1} - a_{j-1} \gamma_{j-1} }
\label{tridiag1}
\\
u'_j &=& \frac{ r_j - a_j u'_{j-1} }{ b_j - a_j \gamma_j } 
\label{tridiag2}
\\
u_j &=& u'_j - \gamma_{j+1} u_{j+1}
\label{tridiag3}
\end{eqnarray}
首先，我们可以证明方程组中的第一个方程被$u_j$满足，即
\begin{equation}
b_1 u_1 + c_1 u_2 = r_1,
\end{equation}
证明如下，根据(\ref{tridiag3})，
\begin{eqnarray}
u_1 &=& u'_1 - \gamma_2 u_2 \nonumber\\
    &=& r_1/b_1 - \frac{c_1}{b_1} u_2,
\end{eqnarray}
即
\begin{equation}
b_1 u_1 + c_1 u_2 = r_1.
\end{equation}

我们可以证明方程组中的最后一个方程，即
\begin{equation}
a_n u_{n-1} + b_n u_n = r_n.
\end{equation}
证明如下
\begin{eqnarray}
a_n u_{n-1} + b_n u_n &=& a_n ( u'_{n-1} - \gamma_n u'_n ) + b_n u'_n \nonumber\\
                    &=& a_n u'_{n-1} + (b_n - a_n \gamma_n)  \nonumber\\
                    &=& r_n,
\end{eqnarray}
上式第一个等号用了(\ref{tridiag3})，第三个等号用了(\ref{tridiag2})。
所以最后一个方程是对的。

对于第$j=2,\cdots,n-1$个方程，将（\ref{tridiag2}）带入（\ref{tridiag3}），得到
\begin{equation}
u_j = \frac{ r_j - a_j u'_{j-1} }{ b_j - a_j \gamma_j } - \gamma_{j+1} u_{j+1},
\end{equation}
而根据（\ref{tridiag3}），
\begin{equation}
u'_{j-1} = u_{j-1} + \gamma_j u_j,
\end{equation}
将上式带入上上式，即得
\begin{equation}
(b_j - a_j \gamma_j) u_j = r_j - a_j (u_{j-1} + \gamma_j u_j) -c_j u_{j+1},
\end{equation}
整理即得
\begin{equation}
a_j u_{j-1} + b_j u_j + c_j u_{j+1} = r_j.
\end{equation}

下面的代码实现追赶法
\begin{lstlisting}
#include<iostream>
using namespace std;

/*
 * TridiagLinear:       Solve linear equations Ax=r, where A has a tridiagonal form:
 *                      b0      c0
 *                      a1      b1      c1
 *                      0       a2      b2      c2
 *                      ...     ...     ...     ...
 *                      ...     ...     ...     ...
 *                                      a_{n-1} b_{n-1}
 * matrix a must be diagonally dominant, otherwise stability problem may appear.
 * The solution is stored in u
 */
void TridiagLinear( int n, double *a, double *b, double *c, double *u, double *r ){

        int j;
        double *gamma = new double [n];
        u[0] = r[0]/b[0];
        gamma[0] = 0;
        for(j=1;j<n;j++){
                gamma[j] = c[j-1]/(b[j-1] - a[j-1]*gamma[j-1]);
                u[j] = ( r[j] - a[j]*u[j-1] )/(b[j] - a[j]*gamma[j]);
        }
        for(j=n-2;j>=0;j--){
                u[j] -= gamma[j+1]*u[j+1];
        }
        delete [] gamma;
}

int main(){
        int n=3;
        double *a = new double [n];
        double *b = new double [n];
        double *c = new double [n];
        double *u = new double [n];
        double *r = new double [n];
        a[0]=0; a[1]=1; a[2]=1;
        b[0]=2; b[1]=2; b[2]=2;
        c[0]=1; c[1]=1; c[2]=0;
        r[0]=3; r[1]=4; r[2]=3;

        TridiagLinear(n, a, b, c, u, r);
        cout<<"u: "; for(int i=0;i<n;i++)cout<<u[i]<<","; cout<<endl;
        
        delete [] a,b,c,u,r;
        return 0;
}
\end{lstlisting}
当系数矩阵$A$为对角占优，即对角元大于非对角元时，追赶法一定会得到方程组的解。

\section{迭代法}
\subsection{雅克比（Jacobi）方法}
线性方程组为
\begin{equation}
A \vec{x} = \vec{b},
\end{equation}
A可以分为三部分，下三角部分（不包括对角元）、对角部分、上三角部分。
\begin{equation}
A = L + D + U,
\end{equation}
所以方程组的解应满足
\begin{equation}
(L+D+U)\vec{x} = \vec{b},
\end{equation}
即
\begin{equation}
\vec{x} = D^{-1} ( \vec{b} - (L+U)\vec{x} ),
\end{equation}
借助这个式子构造迭代式
\begin{equation}
\vec{x}^{(k+1)} = D^{-1} ( \vec{b} - (L+U)\vec{x}^{(k)} ),
\end{equation}
前人证明了，如果矩阵A是正定的，或者对角占优的，这个迭代式一定会收敛至正确解。

\subsection{高斯-席德尔（Gauss-Seidel）方法}
在雅克比方法中，迭代式为
\begin{equation}
\vec{x}^{(k+1)} = D^{-1} ( \vec{b} - (L+U)\vec{x}^{(k)} ),
\end{equation}
这里面实际上有n个迭代式，n为方程组维数。
在等号右边使用$\vec{x}$的矩阵元$x^{(k)}_i$，可能已经产生了第k+1次迭代结果$x^{(k+1)}_i$，那么，如果使用$x^{(k+1)}_i$而不是$x^{(k)}_i$，一般会加速收敛的速度。

这样做，就叫做高斯-席德尔方法。
前人证明了，如果矩阵A是正定的，或者对角占优的，这个迭代式一定会收敛至正确解。

据说，如果不满足这些严格条件，有时Jacobi和Gauss-Seidel方法也都会收敛到正确解。

迭代方法还有逐次松弛法。这里我省略了。。。

\section{应用：最小二乘法}
现在假想有一个物理现象，所有人都很感兴趣，于是实验家、理论家都研究它。

实验家说，最有意思的可观测量是$y$，有$n$个条件都会影响它，$\vec{x}=(x_1,x_2,\cdots,x_n)~$。
所以$y$是$\vec{x}$的函数$y=f(\vec{x})$。
通过几百天夜以继日的实验，我们伟大的实验家测得了$N$组数据，即$(\vec{x}_1,y_1),~(\vec{x}_2,y_2),~\cdots,~(\vec{x}_N,y_N)~$。

理论家说，根据我最聪明的头脑，技巧高超的分析，$y$与$\vec{x}$的函数关系为
\begin{eqnarray}
y = f(\vec{x}) = \sum\limits^m_{i=1} a_i g_i(\vec{x}),
\label{hypothesis}
\end{eqnarray}
其中$g_i(\vec{x})$是已知函数，我可以用手写出来，代表着m个重要的方面。
$a_i(i=1\cdots m)~$是待定参数，可通过拟合实验数据得到。

实验家对理论总是充满怀疑的，他只相信自己的数据。
但他想，就算他们拟合不好，也只会使事情更有趣。

于是他们开始做拟合。
对一组参数$~\vec{a}=(a_1,a_2,\cdots,a_m)$，$y=f(\vec{x})~$在$~N~$个数据点的误差平方和为
\begin{eqnarray}
\Delta(\vec{a})
= \sum\limits^N_{j=1} (f(\vec{x}_j)-y_j)^2
= \sum\limits^N_{j=1} (\sum\limits^m_{i=1} a_i g_i(\vec{x}_j) - y_j)^2.
\label{Delta}
\end{eqnarray}
这是一个关于$a_1, a_2, \cdots, a_m$的二次型，有且仅有一个极小值点。

在极小值处，$\Delta$对任意$a_l$的一阶偏导数都等于零。
即$~\forall l \in 1,\cdots,m~$，
\begin{eqnarray}
\frac{\partial \Delta}{\partial a_l}
&=& \sum\limits^N_{j=1} 2(\sum\limits^m_{i=1} a_i g_i(\vec{x}_j) - y_j) g_l(\vec{x}_j)
\nonumber\\
&=& 2 \left\{ \sum\limits^m_{i=1} (\sum\limits^N_{j=1}g_l(\vec{x}_j) g_i(\vec{x}_j)) a_i - y_j \sum\limits^N_{j=1} g_l(\vec{x}_j) \right\} = 0 .
\end{eqnarray}
即解多元一次线性方程组
\begin{eqnarray}
\left(
\begin{aligned}
~~~~~~~~~~\cdots~~~~~~~~~~~~~~~\\
\cdots~~~
\sum\limits^N_{j=1} g_l(\vec{x}_j) g_m(\vec{x}_j)~~~\cdots\\
~~~~~~~~~~\cdots~~~~~~~~~~~~~~~
\end{aligned}
\right)
\left(
\begin{aligned}
\cdot~~\\
\cdot~~\\
a_m\\
\cdot~~\\
\cdot~~\\
\end{aligned}
\right)
=
\left(
\begin{aligned}
\cdot~~~~~~~~\\
\cdot~~~~~~~~\\
\sum\limits^N_{j=1} g_l(\vec{x}_j)y_j\\
\cdot~~~~~~~~\\
\cdot~~~~~~~~\\
\end{aligned}
\right)
\end{eqnarray}
得到的解$\vec{a}$即最佳参数，将最佳参数带入$\Delta$的表达式，就得到了误差的极小值$\Delta_{min}$。

如果$\Delta_{min}$非常非常小，理论家可能会很得意。
但是如果$m$非常大，即参数个数非常多，实验家可能也不买账，毕竟，给你一百万个参数，大概什么都能拟合好。。。
所以，高贵的理论都是参数少而拟合准的。

下面是实现最小二乘法的c++代码
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>

/*
 * GaussLinear: solves n-dim linear equations
 * The solution is stored in x
 * Permutation of rows is activated, when the principle element is smaller than min
 * If A is singular, a permutation will fail and function exit(1), signals error
 */
void GaussLinear(int n, double **A, double *x, double min){

        int i,j,k;
        double y;
        for(i=0;i<n;i++){
                if( fabs( A[i][i] ) < min ){// permutation is activated
                        for(j=i+1;j<n;j++){
                                if( fabs( A[j][i] ) > min ){
                                        for(k=i;k<n;k++){
                                                y = A[i][k];
                                                A[i][k] = A[j][k];
                                                A[j][k] = y;
                                        }
                                        break;
                                }
                        }
                        if(j==n){// didn't find any element larger than min
                                cerr<<" A is singular in GaussLinear(n,A,x,min).\n";
                                exit(1);
                        }
                }
                for(j=i+1;j<n;j++){// row operations
                        y = - A[j][i] / A[i][i];
                        for(k=0;k<n+1;k++){
                                A[j][k] = y * A[i][k] + A[j][k];
                        }
                }
                // for test only
                /*
                cout<<"i="<<i<<"\t A:"<<endl;
                for(j=0;j<n;j++){
                        for(k=0;k<n+1;k++){
                                cout<<A[j][k]<<",";
                        }
                        cout<<endl;
                }
                */
        }
        for(j=n-1;j>=0;j--){
                y=0;
                for(k=j+1;k<n;k++) y += A[j][k]*x[k];
                x[j] = 1/A[j][j] * ( A[j][n] - y );
        }
}

double LeastSquare(int N, double **x_data, double *y_data, int n, int m, double *a, double (**func)(int, double *), double precision ){

        int l, i, j;
        double y,z;
        double ** A = new double * [ m ];
        for(l=0;l<m;l++){
                A[l] = new double [m+1];
                for(i=0;i<m;i++){
                        y=0;
                        for(j=0;j<N;j++){
                                y += func[l](n, x_data[j]) * func[i](n, x_data[j]);
                        }
                        A[l][i] = y;
                }
                y=0;
                for(j=0;j<N;j++){
                        y += func[l](n, x_data[j]) * y_data[j];
                }
                A[l][i] = y;
        }
        // for test only
        /*
        cout<<"A is ready"<<endl;
        for(i=0;i<m;i++){
                for(j=0;j<m+1;j++){
                        cout<<A[i][j]<<",";
                }
                cout<<endl;
        }
        */
        GaussLinear(m, A, a, precision);
        y=0;
        for(j=0;j<N;j++){
                z=0;
                for(i=0;i<m;i++)
                        z += a[i] * func[i](n, x_data[j]);
                z -= y_data[j];
                y += z*z;
        }
        y = sqrt( y/N );
        for(i=0;i<m;i++) delete [] A[i];
        delete [] A;
        return y;// returns root-mean-square-root
}

double g0(int n, double *x){
        return 1;
}

double g1(int n, double *x){
        double z = 0;
        for(int i=0;i<n;i++) z += x[i];
        return z;
}

double g2(int n, double *x){
        double z = 0;
        for(int i=0;i<n;i++) z += x[i] * x[i];
        return z;
}

double g3(int n, double *x){
        double z = 0;
        for(int i=0;i<n;i++) z += pow(x[i],3);
        return z;
}

int main(){

        int i,j,n=1, m=3,N=10;
        double ( *g[3] )(int, double * )={g0,g1,g2};

        double **x_data = new double * [N];
        for(i=0;i<N;i++) x_data[i] = new double [n];
        double *y_data = new double [N];

        for(i=0;i<N;i++){
                x_data[i][0]=i;
                y_data[i]=1.234+2.456*i+3.789*i*i;
        }

        double *a = new double [m];
        cout<<"RMSD="<<LeastSquare(N, x_data, y_data, n, m, a, g, 0.01 )<<endl;
        cout<<"a: "; for(i=0;i<m;i++)cout<<a[i]<<","; cout<<endl;

        delete [] a, y_data;
        for(i=0;i<N;i++)delete [] x_data[i]; delete [] x_data;

        return 0;
}
\end{lstlisting}
我们定义了$x_{data}, y_{data}$，并人为地设定了10个简单数据点，用于检验新写的代码；
在实战中，你可以用文件流读取数据文件中的数据。
我们假设理论值为
\begin{equation}
f(\vec{x}) = a_0 g_0(\vec{x}) + a_1 g_1(\vec{x}) + a_2 g_2(\vec{x}).
\end{equation}
编译运行上面的代码，你会得到
\begin{lstlisting}
RMSD=1.36473e-13
a: 1.234,2.456,3.789,
\end{lstlisting}
所以得到的最佳参数为$1.234,2.456,3.789$，在最佳参数处得到最小的$RMSD=1.36473e-13$，其中$RMSD$是root-mean-square-deviation，即均方根误差的缩写，
\begin{equation}
RMSD = \sqrt{ \frac{ \sum^N_{i=1} (f(\vec{x}_i) - y_i)^2 }{ N} }.
\end{equation}

在代码中，我们用前文中写的高斯消元法求取线性方程组的根。
如果每次应用，都需要复制粘贴一遍GaussLinear函数，就太费劲了。
我们可以编辑一个头文件，比如叫做library.h，然后把我们定义过的所有函数粘贴在里面，需要用到其中的函数时，就在应用程序中加上一行
\begin{lstlisting}
#include"library.h"
\end{lstlisting}
并保证library.h在同一个文件夹下，即可。